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The self-consistent Gaussian approximation (SCGA) for classical spin systems described by a 
completely anisotropic D-component vector model is proposed, which takes into account fluctuations 
of the molecular field and thus is a next step beyond the molecular field approximation. The SCGA is 
sensitive to the lattice dimension and structure and to the form of spin interactions and yields rather 
accurate values of the field-dependent magnetization m(H, T) and other thermodynamic functions 
in the whole plane (H,T) excluding the vicinity of the critical point (0,T C ), where the SCGA breaks 
down, showing a first-order phase transition. The values of T c themselves can be determined in the 
SCGA with an accuracy better than 1% for actual 3-dimensional structures. At low and high temper- 
atures the SCGA recovers the leading terms of the spin-wave theory, the low- and high-temperature 
series expansions, respectively. The accuracy of the SCGA increases with the increase of the spin 
dimension D, and in the limit D — » oo the exact solution for the spherical model is recovered. 
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I. INTRODUCTION 

The models considering spin as a classical vector vari- 
able of a fixed length are the most studied ones in the 
theory of temperature-induced phase transitions on a lat- 
tice. The reason for this is that quantum effects that are 
always present in real magnetic systems and can make 
calculations more intricate play, however, a secondary 
role at (and above) T c and do not change the critical be- 
havior. Classical models are also a good approximation 
for magnetics with large spin values S, such as, e.g., the 
Heisenberg systems EuO and EuS having S — 7/2. Gen- 
erally, for systems described by the Heisenberg model, 
H(S'), with S > 1 quantum effects become irrelevant in 
the temperature range T > T c /S, where the whole Bril- 
louin zone is populated by spin waves and their occupa- 
tion numbers become large. Additionally, the arbitrary- 
S Ising model, I(S), without a transverse field can be 
treated with the same methods as classical ones, be- 
cause no spin commutators appear and quantum effects 
are trivial. An advantage of classical vector models is 
that they can be formulated for an arbitrMy-jaumber of 
spin components, as was done by StanlejocHa. Such a 
generalization is important since some magnetics with a 
complicated structure possess an order parameter with 
n > 3 symmetric components (see Refs. ^, ^[ and refer- 
ences therein). 

In the absense of an analytical solution to the phase 
transition problem in three dimensions such numerical 
methods as low-temperature series expansions (LTSE's) 
for the 1(1/2) modeE and highjlcm»eratu*e|Series exft&B,- 
sions £H jESE's) for the I(l/2)0BHta, I(5)EM3, H(S)Bii , 
H(oo)li3 : Ej and classical plane rotator and x-y modelsEZI, 
as well as for the general n-component vector model, 



O(n)S'!lll'fi]'0, were successfully applied for an accurate 
calculation of thermodynamic quantities in a wide tem- 
perature range including the vicinity of T c . It gave the 
results for the critical indices of magnetic systems and 
favored the creation of scaling and universality concepts. 
With the development of computational facilities and al- 
gorithms the series methods were permanently improved. 
As the latest benchmark the recent calculationEi of the 
HTSE series for the reduced susceptibility T\{T) of the 
0(n) models up to (J/T) 19 can be considered. Another 
very efficient numerical method competing with series ex- 
pansions is based on Monte Carlo (MC) simulations (see, 
e.g., Refs. ||[ [l(]). An extraction of accurate results for 
infinite systems from the simulation data for the lattices 
with a finite lineac dimension L is based usually on the 
finite-size scalingcJ. An alternative approach also using 
simulations is the chiral perturbation theory in powers of 
1/L (see, e.g., Ref. p4| and references therein). 

Lately the ideas of the statistical theory of magnetism 
together with the methods of calculation have penetrated 
into the field theory. In particular, the lattice-regularized 
scalar Higgs model in the chiral limit, which can be 
identified with the four-componettt-classical Heisenberg 
model, 0(4)w-4n four dimensionsEJt|i|-was studied with 
the HTSEBE3 and MC simtdatiDrSo methods. Re- 
cently Wilczek and RajagopalESEH have related the two- 
flavor quantum chromodynamics (QCD) to the 0(4) vec- 
tor model in three dimensions, the gauge coupling con- 
stant g 2 determining the temperature and the quark mass 
m q being proportional to the applied magnetic field. This 
has initiated an extensive numerical work (see Refs. |2^, 
pc| |30| and references therein) . 

Although HTSE's produce the series coefficients usu- 
ally with the help of such diagram methods as the linked 
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cluster expansion the results are represented 

as a sum of "bare" (unrenormalized) diagrams, each pro- 
portional to some power of J/T. Alternatively, there 
were attempts starting from the early years to sum up 
some "important" infinite diagrams series to obtain a 
closed-form equation for a magnetic system in terms of 
renormalized diagrams, which should be a good analyt- 
ical approximation in the whole temperature range. It 
was shown, in particular, how the mean field approxi- 
mation (MFA) can be obtained diagrammatically (see, 
e.g., Refs. [33], |34|). A further renormalization Cjf-,dia- 
grams for the Ising model by Horwitz and Callent-3 led 
to an improvement of the MFA taking into account self- 
consistently Gaussian fluctuations of the molecular field. 
This important work remained seemingly unappreciated, 
since the resulting equations were not numerically inves- 
tigated in a satisfactory way and the real accuracy of the 
approximation was not recognized. Only much later was 
this self-consistent Gaussian approximation (SCGA) for 
the Ising systems independently rediscovered and numer- 
ically analyzed in Ref. |35|. ■— ■ 

The methods developed by Horwitz and Callcna for 
the Ising model were generalized for the quantum Heisen- 
berg model in the subsequent paper, Ref. |3(| This first 
version of the spin diagram technique (SDT) had not, 
however, succeeded in formulating a Gaussian approxi- 
mation for the Heisenberg systems, since due to quantum 
effects transverse spin cumulants acquire a time depen- 
dence and cannot be renormalized in a desirable way. 

Later a similar diagram technique was formulated in- 
dependently in Refs. [37[ [38] and further developed in Refs. 

Although SDT allows one to write down diagram- 
matic perturbation series for all temperatures in a regular 
way and recovers the known spin- wave and LTSE results 
in the ordered state, as well as the HTSE ones above 
T c , summation of nontrivial diagram sequences in all or- 
ders of a perturbation theory (apart from the usual cases 
of Dyson and vertex equations) seems to be impossible. 
Due to the use of the Wick theorem for the calculation of 
averages of transverse spin components the number and 
complexity of diagrams increase dramatically with each 
order and most of diagrams are divergent at T > T c and 
compensate each other only in final expressions. The lat- 
ter is not the case only for Ising systems, where there are 
no problems with the noncommutativity of different spin 
components and quantum effects are trivial. In the clas- 
sical limit S — > 00, which is of a primary importance in 
the theory of temperature-induced phase transitions, the 
quantum SDT does not essentially simplify. 

In Ref. |4l] an alternative diagram technique for classi- 
cal spin systems was proposed, which explicitly takes ad- 
vantage of their classical properties and is much simpler 
than the quantum SDT. It allows, in particular, calcu- 
lation of thermodynamic quantities of a system without 
dealing with its dynamics. In the static case all spin 
components can be treated similarly, and the consid- 
eration can be carried out for a generalized completely 
anisotropic model of D-component classical spin vectors 
(|m| = 1) on a lattice: 
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Jij T) a m a im a j . (1.1) 



a = l 



If the exchange interaction rj a Jij is isotropic, i.e., all 
anisotropy factors r] a =.l, this model reduces to the one 
pioneered by Stanleyu'ErH, who provedo that it is in the 
limit D -r-* 00 equivalent to the exactly soluble spheri- 
cal mo delEJ . An important particular case of the general 
model (1.1) is the so-called n-D modeled, where n ^ D 



spin components are coupled by the exchange interaction 
with an equal strength and the rest D — n ones are "free" 
(i.e., r} a = 1 for a n and n a — for a > n). The 
n-D model contains as particular cases the S =1/2 Ising 
model, 1(1/2), for n = D = 1, the classical Ising model, 
I(oo), for n = 1, D = 3, the plane rotator model for 
n = D = 2, the classical x-y model for n = 2, D = 3, and 
the classical Heisenberg model, H(oo), for n = D = 3. 
The variable n is the number of the order parameter 
components and determines the universality class of a 
system. The total number of spin components, D, enters 
only such nonuniversal quantities as T c . It is clear that 
the expansion of the critical indices for the large num- 
ber of components can be only the 1/n expansion. To 
the contrary, we shall see below that the absolute values 
of thermodynamic quantities are naturally developed in 
powers of 1/D for D 3> 1, which is not automatically the 
same as 1/n for n ^ D. 

In Ref. |ll| the self-consistent Gaussian approximation 
by Horwitz and Callen was generalized for systems with 
continuous spin symmetry and it was shown that in the 
limit D — * 00 the SCGA becomes exact and yields the 
solution of the spherical model, whereas all other dia- 
grams die out as at least 1/D. Accordingly, the SCGA 
becomes more accurate for high spin dimensions D and 
works better for H(oo) model (n = D = 3) than for 
the 10-/2) one (n= D = 1). Numerical calculations for 
I(5')E3 and H(oo)a models have shown that for differ- 
ent 3-dimensional lattice structures the SCGA yields the 
magnetization m and other thermodynamic quantities in 
the whole temperature range excluding the close vicinity 
of T c with an overall accuracy about 1%, including the 
determination of T c itself. 

In Refs. H, g| H, the SCGA was only briefly de- 
scribed, and its analytical properties need to be explained 
in more detail. Principally important is to test the SCGA 
on models with lattice dimensionality d 4 (hypecubic 
lattices) and_tp compare its results with those of the 1 /d 
expansior£3c3 and MC calculationsEJ. In this case the 
SCGA should be more accurate, since nontrivial effects 



of the 
die out 



ctuation interaction (i.e., non-Gaussian effects) 
ill In view of applications in the field theory men- 
tioned above it is important to extend calculations to 
O(n) models (n = D) with n > A .and to make a com- 
parison with the 1/n expansiorOOCJ. Some other tasks 
are to perform a numerical solution of the SCGA equa- 
tions in the case of a nonzero magnetic field, to make 
a comparison with the experimental data on Eu chalco- 
genides, and to consider the lattices with the next nearest 
neighbor (nnn) interactions. The solution of the prob- 
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lems mentioned above, as well as a detailed statement of 
the SCGA, is the aim of the present article. 

In Sec]lJ a simple derivation and analysis of the SCGA 
for the Ising systems without using diagrams is given. In 
Scc.pl the classical spin diagram technique and construc- 
tion of the SCGA for a general Hamiltonian (1.1) are de- 
scribed in more detail. In Sec. IV the analytic properties 
of the SCGA in different limiting cases are investigated, 
including the spherical limit, where the known results 



are generalized for the anisotropic Hamiltonian (1.1). In 
Sec.[y] the results of the numerical solution of the SCGA 
equations for different classical spin models on different 
lattices are presented and compared with the available 
HTSE, LTSE, MC simulation, and 1/D expansion re- 
sults, as well as with the experimental data on EuO and 
EuS. In Secjv] some further applications of the SCGA 
and the possibilities of its generalization are discussed. 



II. IDEA OF THE SCGA 



which is reflected by the shift of the actual values of T c 
in about 30% downwards from T^ FA depending on the 
l atti ce structure and the details of spin interactions in 
(1.1). The latter makes feasible an improvement of the 
MFA in d ^ 3 dimensions, which consists in taking into 
account molecular field fluctuations described only by the 
set of their second moments o~ 2a . This means that the 
averages of an arbitrary number of molecular field com- 
ponents decay pairwise, which is equivalent to the use of 
the Gaussian distribution function for the molecular field 
fluctuations. For the Ising model (r] a — for a ^ z) this 
leads, in particular, to the expression for magnetization 
m being given by a Langevin function with a spreaded 
argument: 



2a 2z ) 
xB[/3({H z )+H zfi ) 



(2.4) 



If in (1.1) the magnetic field H is directed along the 
ordering axis z (r\ z — 1), then the z component of the 
molecular field IT acting on the spin on a site i is given 
by 



H Z i — H 



(2.1) 



the MFA consists in neglecting fluctuations of IT, which 
in the spatially homogeneous case leads to the Curie- 
Weiss equation for magnetization m = (m z ): 



m = B((3(H z )), 



(H z ) = H + mJ , 



(2.2) 



where B(£) is the Langevin function, f3 = 1/T, and Jo is 
the zero Fourier component of the exchange interaction. 
The second moment of fluctuations of the a component 
of the molecular field IT, which were neglected in the 
MFA , can be expressed as 



<r 2a = 2_j JijJij'vl(Am aj Am af ) 



33' 



= VoJj^ I (v a J< i yS aa (q), (2.3) 

where Am = m — (m z )e z , S aa (q) is the spin-spin corre- 
lation function, vq is the unit cell volume, and d is the 
lattice dimensionality. If correlations of spins on differ- 
ent lattice sites j, f are neglected, then for systems with 
nearest neighbor (equivalent neighb or) interactions the 
integral over the Brillouin zone in ( |2.3| ) is proportional 
to 1/z and small for a large number of equivalent neigh- 
bors z. This is justified in the tempe ratu re range T>T C , 
but for T ~ T c the correlations in (|j^) should be taken 
into account. For low -dim ensional systems (d — 1,2) 
the lattice integral in (2.3) diverges with lowering tem- 
perature at q = 0, which invalidates the MFA . For 
three-dimensional systems the magnitude of the molecu- 
lar field fluctuations a 2 remains finite and not very large, 



m = B{^l z ) 



1-1/2 



dze- z2 B{£ z + 2l 1 z /2 z) 1 (2.5) 



where £ z = (3(H + mJo) and l z = (3 2 a 2z /2. To obtain 
a closed system of equations, one ca n ca lculate the spin- 
spin correlation function S zz (q) in ( ^.3j ) in the simplest 
Ornstcin-Zcrnikc approximation: 



S zz (q) = 



l-B'(£ z ,l z )0Jq 



(2.6) 



where the derivative of the Langevin function, B' = 
dB /d£, is also ren ormalized by Gaussian fluctuations 
analogously to ( |2.5D . This system of nonlinear equations 
for m and l z given by (2.5), ( |2.6| ), and (2.3) with a = z 
was obtained in a very technical manner by Horwitz and 
CallenE.3 and was solved numerically in Ref. p35|. Note 



that the integral over t he B rillouin zone a 2z , Eq. (2.3), is 
taken into account in (2.5) in all orders of a perturbation 
theory. Such a self-consistent Gaussian approximation 
is, like all closed-form approximations in the theory of 
phase transitions, not a rigorous expansion in some small 
parameter. It is an approach taking into account some 
physically significant diagram structures self-consistently 
in all orders of a perturbation theory and reproducing 
the leading orders in the perturbatively treatable regions 
r«T c and T > T c . In the next section the SCGA will 
be derived for a general form of the spin-vector Hamil- 
tonian (1.1) with the use of the classical spin diagram 
technique. 



III. CLASSICAL SPIN DIAGRAM TECHNIQUE 
AND THE SCGA 

This diagram technique can be considered as a siriu 
plified form of the quantum linked cluster expansionL3 
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or of the quantum SDT0S, making use of the classical 
properties of spin vectors. A perturbative expansion of 
the thermal average of any quantity A characterizing a 
classical spin system (e.g., A = m z ) can be obtained by 
rewriting ([0]) as H = Ho + H m t , where Ho is the MFA 
Hamiltonian with the molecular field (H z ) determined by 
( |2.2| ), and expanding the expression 

i r N 

3=1 



in powers of 7ii n t ■ The integration in ( 3T ) is carried out 
with respect to the orientations of the _D-dimensional unit 
vectors nij on each of the total N lattice sites. Averages 
of various spin vector components on various lattice sites 
with the Hamiltonian Ho can be expressed through spin 
cumulants, or semi-invariants, which will be considered 
below, in the following way: 

{m ai )o = A„, 

(m ai m/3j}o = k a p8ij + A a Ap, (3.2) 
(m ,impjm 1 k)o = Aa^Sijk + A a pAj6ij 

+ AfjjA a 5j k + AyaApSki + A ct A / 3A 7 , 

etc., where <5y, Sijk, etc., are the site Kronecker symbols 
equal to 1 for all site indices coinciding with each other 
and to zero in all othe r cases. For the one-site averages 
(i = j = k = . . .) (3.2) reduces to the well-known repre- 
sentation of moments through semi-invariants, general- 
ized for a multiple-component case. In the graphical lan- 
guage (see Fig. l|) the decomposition (3.2) corresponds to 
all possible groupings of small circles (spin components) 
into oval blocks (cumulant averages) . The circles coming 
from Hi n t (the "inner" circles) are connected pairwise 
by the wavy interaction lines representing the quantity 
r] a pjij in (1.1). In diagram expressions summations over 
site indices i and component indices a of inner circles 
are carried out. One should not take into account dis- 
connected (unlinked) diagrams [i.e., those containing dis- 
connected parts with no "outer" circles belonging to A in 
(3.1)], since these diagrams are compensated for by the 
ex pans ion of the partition function Z in the denominator 
of (3.1). Consideration of combinatorial numbers shows 
that each diagram contains the factor l/n s , where n s is 
the number of symmetry group elements of a diagram. 
The symmetry operations do not concern outer circles, 
which serve as a distinguishable "root" to build up more 
complicated or renormalized diagrams. Such combina- 
torial factor s are present, in particular, in the formulas 
( 3.15 ) and ( 3.16 ). For practical calculations it is usu- 
ally more convenient to use the Fourier representation 
and to calculate integrals over the Brillouin zone rather 
than lattice sums. As due to the Kronecker symbols in 



(3.2) lattice sums are subject to the constraint that the 
coordinates of the circles belonging to the same block 
coincide with each other, in the Fourier representation 
the sum of wave vectors coming to or going out of any 
block along inte raction lines is zero. The cumulant spin 
averages in (|3.2j ) can be obtained by differentiating the 



generating function A(£) over appropriate components of 
the dimensionless field £ = /3HEJ: 



A aia2 ... ap (£) — 



0"A(f) 



A(0 = lnZ (0, 
where £ = |£|, 

2o(0 = const xr (D/2_1) lD/2-i(0 



(3.3) 



(3.4) 



is the partition function of a /^-component classical spin, 
and I„(£) is the modified Bessel function..—. A similar 
technique was applied by Lusher and Weiszo to gener- 
ate HTSE's for a more general n-component </> 4 m odel . 
For several lowest-order cumulants differentiation in ( |3.3| ) 
leads to the following expressions: 

Aa(*) = Bo(0&»=-B(OW6 

A Q/3 (0 = B (0 S a /3 + £a£/3, (3.5) 

A Q /3 7 (£) = (£a Vr + £/3<5 7 a + £~/S a f3) 

Aap-ysiO = Bi 3V(S a pS^s) 

+ B 2 6V(^6 y5 ) + B 3 te&ts, 

where 8 a p is the spin component Kronecker symbol, V is 
the symmetrization operator, 



ld_\ n B(£) 



and 



B(0 = dA(0/d£ = I D/2 (Ofin/ 2 -i(0 



(3.6) 



(3.7) 



is the Langevin function of D-component classical spins, 
which can be expressed through elementary functions for 
odd values of D: 

tanh(£), D = 1, 

B(£) = < coth(0 - D = 3, (3.8) 

l/(coth(0 - 1/0 - 3/£, D = 5, 

etc. The small- and large- argument expansions of the 
Langevin function £?(£) have the form 



B(0 



e 



2£ 5 



D D 2 {D + 2) D 3 (D + 2)(D + 4) 
5£ 7 ( ^ 2 



(3.9) 



D 4 (D + 2)(D + 4)(D + 6) 



1 



5(D + 2) 



and 



B(0 S 1 - 



D-l (D-l)(D -3) 



2£ 



(3.10) 



respectively. On e can see from ( p.9[ ), that the functions 
B n (£), Eq. Q, are all finite at £ = 0: B (0) = 1/D, 
B^O) = -2/[D^(D + 2)], B 2 = 16/[£» 3 (D + 2)(£> + 4)], 
etc. Accordingly the spin cumulants A... in (3.5) with an 
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even number of coinsiding indices are given in this case 
by their first terms: 

A Qa = Bq(0), A aa p =B 1 (O)(l + 2S a0 ), (3.11) 



etc., whereas all other cumulants 



turn to zero. At 
follows B n (£) oc 
5|) yield compara- 



large arguments from (3.6) and ( |3.1 
^-(i+2n)_ In this limit all terms f(| 

ble contributions into A..., and a fc-spin cumulant decays 
generally as A Q1Q2 ... Qfc oc £ _ ( fc_1 ). If, however, the field £ 
is directed along some axis z, then in the cumulant aver- 
ages containing z components of spins the leading terms 
can cancel e ach other. In particular, the two-spin cumu- 
lant A a p in ( |3.5| ), which plays a big role in the following, 
can be rewritten explicitly as 



A Q/3 (0 



S a a — 



£a£/3 

e 



B'(0 



(3-12) 



For £ = £e z this expression simplifies to A 22 = -B'(£) 
and A aa = £?(£)/£ (a ^ z). Now from ( 3.10 ) one can see 
that, for £ 3> 1, A 22 cx £~ 2 , whereas A aa oc £ _1 . 

The simplification of spin cumulants for £ = £e 2 men- 
tioned above takes place in the unrenor mali zed diagrams 
generated initially by the expansion of (3.1) in powers of 
7ii nt since there is only one nonzero component of the 
molecular field: £ 2 = £ = (3{ H + m Jp). The complete 
form of spin cumulants (3.E), ( |3.12] ) is needed, however, 
for the construction of the SCGA, which allows for both 
longitudinal and transverse fluctuations of the molecular 
field. The latter is the essense of the diagram technique 
for the multiple-component classical spin systems pre- 
sented here. In the Ising case the classical spin diagram 
technique coincides -with—the "Ising part" of the stan- 
dard quantum SDTt3'E§E3 and can be used with Bril- 
louin functions Bg of a general spin S. In the Refs. |3^, 
|50| the reader can find more technical details concerning 
the construction of SDT for Ising systems, which play 
the same role in the present classical SDT. 

Before proceeding to the construction of the SCGA we 
should make a remark about the numerical calculation of 



the generalized Langevi n fu nction B(£) (3.7) for arbitrary 
D. One can see from ( j3~§| ) that for D > 1 the function 
-B(£) contains terms divergent at £ — > 0, although -B(£) 
itself is well behaved. This hampers numerical calcula- 
tions, and the situation is aggrava ted for the derivative 
£?'(£) and f or th e functions -B n (£) ( |3.6| ) entering the spin 
cumulants (3.5), as well as for higher spin dimensional- 
ities D. The best way of calculating -B(£) is based on 
using the backward recursion relation with respect to D: 



B(D,0 = 



£ 



D + ££(£> + 2, £) : 



(3.13) 



which can be derived from (3.7) and the three-term re- 



cursion relation for the modified Bessel functions Ii/(£)- 
This formula yie lds t he proper small-argument behavior 
of B(D,^), Eq. (|3.9|), to leading order even for an inac- 
curate B(D + 2, £), and the prop er behavior at £ ;§> 1 to 
leading order described by ( 3.10| ) can be guaranteed, -if 
we choose the first two terms of the large- 1? expansion^ 



t JVc^ 



AAA/V = A/WV + 



WW 




+... 



FIG. 1. Self-consistent Gaussian approximation (SCGA) 
for classical spin systems, (a) and (c): block summations for 
the renormalized magnetization and pair spin cumulant av- 
erages; (b): Dyson equation for the renormalized interaction 
line. 



B(o^m + ^f( X ) + o^ 

x = 2£/A f(x) 



1 + (1 + x 2 ) 1 / 2 : 



(3.14) 



as an initial condition for the recurrence formula ( 3.13[ ). 
This procedure pro ves t o be extremely good: Already 
one application of ( 3.13| ) yields -B(£) with an accuracy 
not worse than 0.6% for D = 1, 0.35% for D = 2, and 
0.25% for D — 3 in the whole range of £, and the process 
converges fast with the increase of the iterations number. 

The self-consistent Gaussian approximation consists 
in taking into account pair correlations of the molec- 
ular field acting on a given spin from its neighbors, 
which implies the Gaussian statistics of molecular field 
fluctuations. The corresponding diagram sequence is 
represented in Fig. [I] and is equivalent to the follow- 
ing closed system of nonlinear equations for magneti- 
zation m = (m z ) and the n orm alized second moments 
l a = /3 2 ct 2q /2! [cf. Q and (|J)] of the molecular field 
fluctuations: 



m = A 2 , 

1 f dq VaPJq 

In, = -TTV 



(3.15) 

2r7(2^1-A QQW 3J q ' a = 1 ^--°- 

Here the spin cumulants A... renormalized by Gaussian 
fluctuations of the molecular field are given according to 
Figs. [J(a), |l|(c) by the series 



A... = A... + ^A... QQ Z Q (3.16) 

a=l 



a, (3=1 



whe re t he "bare" spin cumulants A... are given by (3.3) 
or (3.5). This series describing the influence of pair- 
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correlated fluctuations of different components of the 
molecular field can be rewritten as 



A 



a=l n a =0 
D 

exp 



Ti 



A...(€). 



(3.17) 



Such exponential differential operators were considered 
by Horwitz and Callena for the Ising model. A gener- 
alization of their results for the multiple-component case 
yields the closed formula 



A = 



1 



,D/2 



"A...(C), 



where £ is the spreaded molecular field given by 



C = (3(H + rnJ Q )e z 



a=l 



(3.18) 



(3.19) 



e Q is the unit vector in the direction a, and the inte- 



gration in (3.18) is performed with respect to the D- 
component vector variable r = {r a }. It can be seen 
that the renormalized spin cumulants A are functions 
of m and all l a . In the Ising case the SCGA system of 
equations fcl5| ) reduces to the one obtained by Horwitz 
and Callena, which was de scribe d in the preceding sec- 
tion, since here only l z in ( [3.15 ) is nonzero and K zz = 
B'(m , l z ). The expression for l a in (3.15) differs from 
(2.3) in that a zero term of the type J dqJ q ~ Ju = 



was added for convenience, which allows one to formulate 
the diagram technique in terms of renormalized interac- 
tions. 

The number of unknown va riable s in the nonlinear sys- 
tem of the SCGA e qua tions ( 3.15 ) is for a general form 
of the Hamiltonian (1.1) equal to D + l. Thus, for exam- 
ple, for a completely anisotropic Heisenberg model there 
are four unknown variables: m, l x , l y and l z . In a more 
complicated case with the magnetic field transverse to 
the ordering axis z, which is not considered here, one 
should take into account different magnetization compo- 
nents and nondiagonal moments of molecular field fluc- 
tuations, lap with a 7^ /?, and the number of unknowns 
increases considerably. Similar takes place for antifer- 
romagnets in magnetic field and for multiple-sublattice 
structures. Of a practical interest is the case when in the 
Hamiltonian (1.1) there are groups of equivalent trans- 
verse (a ^ z) spin components having anisotropy factors 
rj a , and hence the moments l a , equal to ea ch other. In 
this case the number of unknowns in ( 3.1 5| ) diminishes; 
denoting such a group with the index x and introducing 



and 

Tlx 

the number of equivalent com- 



ponents in t he gr oup, one can simplify the D-dimensional 
integral in (3.18) with the help of the identity 



_J_/"d- 



T(n x /2) 



dr x r™* V 



(3.20) 



and make a replacement =>• ^/u x , where £ x = 2l x ^ 2 r x 
and l x = l a exi in t he pair spin cumulants A QQ , Eq. 



(i 



12), entering ( 3.18 ). Thus, in particular, for the 0(n) 
model all components with a ^ z are equivalent, and 
there are th ree independent variables in the SCGA equa- 
tions ( |3.15| ): m, l z and l x = l a . The Gaussian integrals 
(fr|) reduce in this case to two-dimensional ones over r z 
and r x , and in ( |3.20 ) n x = n— 1. Above T c in the absence 
of a magnetic field m = and all spin components are 
equ ivale nt; there is only one unkno wn va riable l z = l a 
in ( [3.15 ), and the intergal A aa , Eq. ( 3.18 ), becomes one 
dimensional. 

The SCGA system of equations ( 3.15| ) determines the 
equation of state of a magnetic system, i.e., the magne- 
tization as a function of temperature and magnetic field, 
m(T,H). The caloric properties of a magnetic system 
in the SCGA can be also determined. In particular, the 
energy of a spin system U = C H) can be obtained by av- 



eraging the Hamiltonian (1.1) and using the expression 
for the renormalized spin correlation function S aa (q) de- 
termined by (2.3) in the form 



S'aa(q) = 



1 A-aaVotft Jq 



(3.21) 



[cf. ( |2.q )1. Using the definition of l a in ( 3.15 ) one gets 

U = -Em - ^J m 2 - ^oJj^jdJq ^ VaS aa (cj) 

1 ° 
= -Hm- -J m 2 -T^2l a A aa , (3.22) 

a=l 

i.e., the energy can be obtained as a by-pr oduct of the 
numerical solution of the SCGA equations ( 3.15| ). Now 
the heat capacity Ch = dU (T, H) jdT can be obtained 
by the differention of fl3.22| ). The most strong result 
is, however, that for the free energy F — ~T\nZ of a 
system. Its diagrammatic derpation, which was accom- 
plished by Horwitz and CallerO for the Ising model, is a 
rather complicated combinatorial problem, since the free- 
energy diagrams have no distinguishable outer circles, 
which could be used as a root for building renormalized 
diagrams. But the generalization of the corresponding re- 
sults for the multiple-component case is straightforward 
and yields 



(3F = ^ J m 



A 



D 

E 



D 

^ ^ la Aaa ; 

a=l 



(3.23) 



where A is the generating function of spin cumulants ( |3.3| ) 
renormalized by Gaussian fluctuations [see ( 3.18 )] and 



(3.24) 



Considering in ( |3.23| ) m and l a as free parameters, 
i.e., F — F(T, H,m, {l a }), and using the identities 



G 



<9A.../<9£ Q = A... a and dA,,Jdl a = A,,, aa , one can ob- 



tain the SCGA system of equations ( 3.15 ) from the re- 
quirement that F be stationary with respect to m and 
I a- dF/dm = and dF/dl a = 0. The expression f or th e 
energy U, Eq. ( 3.22 ), can be also obtained from ( 3.23 ): 
U = d([3F)/d(3. 



IV. ANALYTICAL PROPERTIES OF THE SCGA 
AND THE SPHERICAL LIMIT 

In this section the behavior of the SCGA solution for 
classical spin systems is analyzed in the regions of high 
and low temperatures, in the spherical limit (D — > oo) 
and in the vicinity of the critical point. It is convenient 
to choose the dimensionless temperature variable 9 = 
T/T C MFA , where t mfa _ j q j d ^ and the dimensionless 

magnetic field h = H/Jq, and suscepti bility x = JoX- 
Then the (unspreaded) molecular field in ( |3.1£ ) is written 
as £ z — f 3(H + m Jo) = (D/8)(h + m), and the quantities 
l a , Eq. (3.15), transform to 



In — ~ [P{f]aGa') 1], Ga — Q A aa , 



20G a 



dq 1 



(4.1) 



where A q = Jq/Jo satisfies 1 — A q oc k 2 for a$k <C 1; ao 
is the lattice spacing. The lattice integral P(X) has the 
following properties: 

l + A 2 /z, A<1, 
P{X)^ { W -c{8Xfl 2 , SX<l,d = 3, (4.2) 

W - cSX\n(c'/SX), <5A<1, d = 4, 



where (5 A = 1 — A, z is the number of equivalent 
neighbors and W (the Watson integral) and c, c' are 
lattice-dependent constants. For low-dimensional sys- 
tems (d = 1, 2) the function P(X) diverges for X — ► 1; for 
<i ^ 5 the leading term of the expansion of P(X) about 
X = 1 is nonsingular. The values of the Watson integral 
W are 1.34466 for the fee lattice (z = 12), 1.39320 for 
the bec lattice (z = 8), 1.51639 for the sc lattice (z = 6), 
1.79288 for the diamond lattice (z = 4), 1.23965 for the 
d = 4 hypercubic (hpc) lattice (z = 8), and 1.15631 for 
the d — 5 hpc one (z = 10). For hypecubic lattices with 
d » 1 one has VF = 1 + 1/z with z = 2d. The differ- 
ence W — 1 measures in the SCGA deviations from the 
molecular field behavior and tends to zero if z — > oo. 

In the high-temperature region (9 3> 1) the second 
mom ents of the molecular field fluctuations, <Jia, Eq. 
(2.3), should be temperature-independent and, corre- 
spondingly, l a = 2 <J2a/2\ oc #~ 2 . In this case the 
renormalized spin cumulants A... ( p. 18 ) are given for 
l < 1 by the expansion J3.16] ), where the 1-st order 
terms written out correspond to diagrams wit h on e in- 
tegration loop in Figs. |l|(a), |l|(c). Now using (3.11) one 
can calculate the quantity G a in fl4.1| ) in the lowest or- 
der: G a S 0- 1 < 1. Then from Q) and (gj) follows 



^« — r i a P > /{^ 2z ) ^ 1; which justifies the initial assump- 
tion. The latter can be used to find a more accurate value 
of G a up to 8~ 3 from the first two terms of (3.16) and 
( |3.1l|) . This allows one to determ ine th e reduced suscep- 



tibilities 9x a {q) = DS aa {q) [see fl3.2l| )] in the SCGA up 
to 8~ 3 . For a particular case of the n-D model (i] a = 1 
for a ^ n and r\ a — for a > n) we write down the 
complete expression for the longitudinal (a n) reduced 
susceptibility up to 8~ 3 , which can be obtained with the 
help of the classical SDT without using the SCGA and 
has the form 



1 

1 + 7T 



_ l n + 2 \ _L_ 

z D + 2 ) 8 2 



2 n 
1 



zD 



1 



(4.3) 



Here all terms except the last one are contained in the 
SCGA, the latter being relatively small as l/[z(D + 2)]. 
Such a situation takes place in the high-temperature 
range for other thermodynamic functions [e.g., the energy 
Q3.22D ], too, as well as in higher orders of a perturbation 
theory - corrections to the SCGA are determined by two 
sma ll parameters 1/z and l/(D + 2). It can be seen from 
(4.3) that for the models with n = D the dependence 
of x\\ 011 D comes practically only from these correction 
terms and remains weak not too close to T c . In the SCGA 
the D dependence of the reduced susceptibility of a spin 
system, as well as of its energy, appears only in the order 
9~ 7 due to the last 1/(D + 2)-correction term in ( |3.9| ) 
and is very weak. For this reason also the values of T c 
determined in the SCGA from the divergence of suscep- 
tibility are for the models with n = D very close to each 
other and to the one of the spherical model. From the 
expression ([0]) it can be seen that in the case of a large 
number of spin components the susceptibility developes 
in a natural way in powers of 1/D and not of 1/n, as 
was mentioned in the Introduction. The same is valid for 
other thermodynamic quantities as well. 

In the low-temperature region (9 C 1) the expansion 
of thermodynamic quantities in powers of 9 is more com- 
plicated, because the longitudinal and transverse spin 
components are nonequivalent and all expressions de- 
pend on magnetization, which should be calculated self- 
consi stent ly in each order. The small-fluctuation expan- 
sion ( 3.16| ) is valid in the range 6 <C 1, too, since the 
high-order spin cumulants diminish as ap propr iate pow- 
ers of l/£ oc 8 [sec the discussion after (3.11)]. In the 
zero-field case, starting from £ = (D/9)m = D/9, one 
can esti mate different terms of the low-temperature ex- 
pansion (frl6b for the magnetization m = A 2 . One gets 
(a ^ z) 

K Z =B = 1-{D- l)/(20 £ 1 - 9{D - 1)/(2D), 
Acta = B/i = 9 jD, 



A 2 
A, 



B' = (D — l)/(2£ 2 ) ^ 9 2 {D - 1)/(2D 2 ), (4.4) 
= (d/dO(B/0 = -9 2 /D 2 , 



A zzz = 9 3 (D - 1)/D 3 , 
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etc., the first of these formulas being the MFA expres- 
sion for magnetization m up to first order in 9. Now 
it can be seen t hat in the low-temperature range G z = 
(D/9)A ZZ oc 9, (4.1) yields l z ~ const, and the contribu- 
tion _of longitudinal fluctuations into m given by A zzz l z 
in ( 3.16 ) is small as 9 3 . The leading contribution to m 
comes from transverse fluctuations (spin waves), since 
G a = 1 and l a = [P{rj a ) - 1] D/{26) > 1. For the ab- 
solute value of the second momen t of the molecular field 
fluctuations, a 2a — 2T 2 l a , Eq. (2.3), the latter means 
02a oc S C 1. Now the magnetization m is given for 
9 <C 1 by the formula 



D 



a=2 



(4.5) 



of the lowest-order spin-wave theory, where the sum in- 
cludes only transverse components. For the n-D model 
(4.5) simplifies tccil 



mea-_[(„ 



l)W + D-n] 



(4.6) 



[see (4.2)]. Such a linear dependence replaces for classical 
fcrromagnets the quantum Bloch law m = 1 — a9 3 ^ 2 . 
In the next order of a perturbation theory in 9 <C 1 



with the help of (3.16) and (4.5) one gets 



Ga = 1 — (6/D)[P(rj a ) — 1], 



(4.7) 



which in the case rj a — 1 leads for three-dimensional sys- 
tems due to (4.2) to the singular negative contribution to 
l a , Eq. (hi), and, as a consequence, to a positive contri- 
bution to magnetization oc 9 3 ^ 2 in addition to the leading 
negative linear term in (4.5). The latter is an artifact of 
the SCGA related to the unbalanced renormalization of 
spin-spin correlation functions. This is, however, an ef- 
fect of the next order of magnitude, which is suppressed 
by the magnetic field or in the anisotropic case r\ a < 1. 

In the spherical limit D — > oo, the SCGA becomes 
exact,— since all other more complicated diagrams die 
outcMil as at least l/D. The Langevin function (3.7) 



simplifies in this limit to the first term of the formula 
( jOg ). The expression for the square of the spreaded 
value of the scaled argument x = 2(/D in ( |3.18] ) reads 



16 



D 

a=2 



(4.8) 



Since, according to (4.1), l a , z oc D, the spreading of the 



z component of the molecular field in (4.8) can be ne- 



glected, whereas the transverse contributions to (4.8) 



each of them is small as l/D, are essential due to their 
number of the order D. The renormalized cumu lants 
A QQ (a ^ z) entering the SCGA equations ( 3.15| ) are 
in the limit D — > oo all equal to each other and given 
according to (Rja) by A QQ = B , so that we can intro- 
duce G = (D/9)Bq. The Gaussian integrals (3.18) are 



for D > 1 easily calculated by applying the identity 



1-1/2 



dx e~ x f{ax 2 ) /(a/2), a < 1, (4.9) 



for an arbitrary function / successively D — 1 times. 
Thus, the integration leads simply to the replacement 
r 2 a => 1/2 in Q, and the SCGA equations ( |U5| ) re- 



duce after some transformations to 



m 

T 



G 



and 



m 



l-G 



D 



a=2 



(4.10) 



(4.11) 



Comparing these results with (4.5) one can identify the 
spherical model as a model which is described in the 
whole temperature range by an effective lowest-order 
spin-wave theory. In a zero magnetic field (h = 0) the 
magnetization m disappears abov e T c , and the quantity 
G, which can be determined from (4.11), increases from 
to 1 with lowering temperature from oo to T c . Below T c 
for h = fr om ( 4.1 0| ) follows G = 1, and m 2 determined 
from (4.11) is a linear function of temperature. It turns 
to zero at the critical point 



q(«0 — 



(4.12) 



which reduces to the well-known resultH 6*e°°' 1 = 1/ 
in the isotropic case r\ a = 1 considered by Stanle; 
The corresponding result for the n-D model (n/D = 
const^ 1) was obtained in Rcf. Wtt. The gener al for 



mula ( 4.12 ), as well as the whole equation of state ( 4.10 ), 
(4.11), shows a crossover to the MFA behavior in the case 
rj a — > 0, rjz = 1 , i.e., for the "spherical Ising model"; see 
(1.1). In the anisotropic case, i.e. for any r\ a < 1, the 



singularity of the function P(X) at X = 1 [see (4.2)] is 
suppressed, and the critical indices of the spherical model 
coincide with those of the MFA . 

Now we proceed to the investigation of the behavior of 
the SCGA solution for classical spin systems in the crit- 
ical region. The first step is to search for T c as a point 
at wh ich the longitudinal correlation function given by 
(|3.21 ) with a = z diverges at q = for h = 0. This 
leads to the condition G z = (D/9)A ZZ = 1, which in the 
isotr opic c ase (rj a = 1) with the use of the symmetriza- 
tion (3.20) canJae transformed to the following nonlinear 
equation for 9JE1: 



2D 



T(D/2) 



dr r D - x e- r F(2l 1 J 2 r), 



l c = D(W - l)/(20 c ), 



where [cf. (|3.12|) 1 



HO = 1 



1 \ B(0 , B'(0 



D 



D 



(4.13) 



(4.14) 
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In the particular case D — 1 this equation reduces to 
the one obtained by Horwitz and CallenO for the Ising 
model. As was stressed above by the analysis of the sus- 
ceptibility HTSE, Eq. (4.3), this 9 C should be very close 
to that of the spherical model, the latter underestimat- 
i ng T c in (5-8)% for three-dimensional systems. Equation 
(4.13) can be solved analytically in two limiting cases: (i) 
for D 3> 1 using the 1/D expansion results of Ref. 51 and 
(ii) for W — 1 « C 1, w hen the spreading of molecular field 
fluctuations in ( 4.13|) is small and the deviation fro m th e 
spherical result is due to the last correction term in (3.9). 
In these limiting cases 9 C is given by 



1 - 



W 



2 (W- l) 3 
D W(2W-l) 2 ' 



D > 1, 



D + 2 



(W-l) 3 , W-Kl. 



(4.15) 



Consi deri ng the values of the Watson integrals W listed 
aft er (|4.2|) , one can see that, indeed, the correction terms 
in ( 4.15] ) are typically small. 

An attempt to simplify the SCGA equations ( 3.15| ) 
about such a defined transition temperature 9 C and to 
calculate the spontaneous magnetization m just below 
9 C shows that 9 C is actually the lower spinodal boundary 
of a fictitious first-order phase transition occurring in the 
SCGA due to its inaccuracy in a close critical region; i.e., 
the magnetization jumps to a finite value by crossing 9 C 
from above. This instability is due to the si ngu lar be- 
havior of the function P(X) near X = 1 [see (4.2)]: The 
decrease of G z from 1 below 9 C related to the increase of 
magnetization leads to a sharp decrease of molecular field 
fluctuations and hence to a further increase of magnetiza- 
tion and so on. Analytically the absence of a continuous 
solution m(9) below 9 C can be shown the most easily for 
the Ising model, where in the vicinity of 9 C the SCGA 
simplifies to a system of equations 



6G Z 
5G Z 



(2D/9)B"' Sl z = 2[(D/8)& - 1], 

{\/{D + 2))(D/9fB"'m 2 = 0. (4.16) 



Here the spreaded derivatives B'") are calculated with 
Lc = D(W — l)/(20 c ) and Sl z is determined as 5l z = 
he - D[P(1 - 5G Z ) - 1]/[20(1 - 5G Z )\ > 0. Below 9 C 
the right part of the first of Eqs. (4.16|) is positive, and 



this equation has no solution since the negative singular 
term with 6l z (B'" < 0) dominates over the positive one 
with SG Z . This is the case for lattice dimensions d = 
3,4; for d ^ 5 the situat ion depends on the numerical 
factors in (4.16), and if ( 1.16| ) has a solution, then the 
MFA behavior of the spontaneous magnetization with 
the critical index (3 — 1/2 is reproduced. The latter_is 
consistent with the analysis of Larkin and Khmelnitskic3. 

The breakdown of the SCGA in the close critical region 
shows that this approximation is more sensitive to criti- 
cal effects than other closed-form approximations always 
reproducing the MFA behavior. In the next section it will 
be shown that the upper spinodal boundary of the phase 
transition determined from the temperature dependence 



of magnetization yields a much better approximation for 
T c than the lower one. It is so because the inverse sus- 
ceptibility turns to zero at T c with zero derivative, and 
even small inaccuracies in determination of the suscep- 
tibility can exert a great effect on determination of T c . 
On the contrary, inaccuracies in magnetization produce 
a smaller effect on T c due to the infinite slope of m at T c . 



V. NUMERICAL RESULTS AND COMPARISONS 



The SCGA system of nonlinear equations (3.15) was 
solved for differe nt l attices and different types of the 
spin hamiltonan (1.1) by the Newton- Raphson iterative 
method. For the fee, bec, sc and diamond lattices the 
analytic expressions for the lattice integrals P(X), Eq. 
(4.1), given by Joyceo were used. For hypercubic lat- 
tices P(X) was reduced to one-dimensional integrals with 
the modified Bessel function Io(a;) and calculated nu- 
merically. In other cases P(X) was calculated by a di- 
rect integration over the Brillouin zone (k XtVyZ G [0,7r]) 
with the use of three-dimensional product quadratures 
composed of 5- or 10-point one-dimensional Gaussian 
quadratures. The accuracy of these quadratures is so 
high that one does not need to a naly tically separate 
the divergence of the integrand in (4.1) at q — » for 
X = 1. The Gaussian integrals ( 3.18] ) were calculated 
for the Ising model with the use of the 5-, 6- or 8-point 
Gauss-Hcrmitc quadratures, for the x-y, plane rotator, 
and completely anisotropic Heisenberg models, with the 
use of the corresponding product quadratures. For the 
models with equivalent spin components, such as O(n ) 
with n ^ 3, the symmetrized integrals of the type ( 3.20| ) 
were approximated by the 5-, 6- or 8-point generalized 
Gauss-Hcrmitc quadratures corresponding to the weight 



function |x| a exp(— x ) with a = 1,2,3 (Ref. |53J) and 
a = 2,4,6,8 (Ref. [34j). The latter was sufficient to cal- 
culate O(n) models up to n = 10. The relative accuracy 
of calculations is not worse than 0.1%, which exceeds the 
intrinsic accuracy of the SCGA. 

The results represented in Table I, Fig. || and Fig. || 
show that the deviations from the MFA due to molecular 
field fluctuations increase with the inverse of the number 
of interacting n eigh bors, z, or, rather, with the differ- 
ence W — 1 [see (4.2)] depending on the lattice structure. 
Among the three-dimensional lattices considered here the 
extreme cases are the diamond lattice (z = 4) and the 
equivalent neighbor fec-sc lattice (z — 18) with 12 face- 
centered-cubic nearest neighbors and 6 simple cubic next 
nearest ones. For the O(n) models (n — D) the devi- 
ations of the magnetization m from the MFA solutions 
increase with the increase of n (see Fig. ||): 1(1/2) => 
plane rotator H(oo) spherical model, whereas the 
susceptibilities x are practically the same for all models. 
The latter could be expect ed from the analysis of the 
susceptibility HTSE, Eq. (4.3), and of the l ower spin- 
odal boundary of the SCGA equations (4.13). For the 
n-D models with increasing n and D — const [I(oo) => 
classical x-y => H(oo)] the deviations from the MFA are 



9 




SCGA for O(n) models 
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FIG. 2. Temperature dependencies of spontaneous magne- 
tization and zero-field susceptibility of the O(n) models on 
the sc lattice in the SCGA. 
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FIG. 3. Temperature dependences of spontaneous magneti- 
zation and zero- field susceptibility of the S — 1/2 Ising model 
on d-dimensional hypercubic lattices in the SCGA, compared 
with Pade-approximations of Refs. rcA H for d — 3. 



increasing stronger than for O(n) ones. This feature is in 
accordance wi th t he functional form of the susceptibility 
HTSE's, Eq. Q . 

The temperature dependences of magnetization and 
other thermodynamic quantities calculated with increas- 
ing temperature are smooth functions of T up to some 
"upper spinodal boundary" after which in a zero mag- 
netic field magnetization jumps to zero. This feature 
results from the inaccuracy of the SCGA in a close criti- 
cal region and was discussed in more detail at the end of 
Sec. IV. The discontinuities of thermodynamic functions 



in the SCGA diminish with the increase of the number of 
interacting neighbors, z, and the number of spin compo- 
nents, D, as well as with the decrease of the number of 
interacting spin components, n. For the quantities which 
are less singular at the critical point (e.g., the energy; see 
Ref. [l3]) these discontinuities are essentially weaker than 
the corresponding deviations from the MFA described by 
the SCGA. 

As was mentioned at the end of the preceding section, 



the SCGA upper spinodal boundary should provide a 
good estimate of phase transition temperatures of three- 
dimensional systems. Indeed, the corresponding values of 
6 C = T C /T C MFA listed in Table I differ from the ones ob- 
tained by HTSE and other accurate methods generally by 
(1-0.5)%, which was never achieved by some other closed- 
form approximation. One can see that for the models 
with higher spin dimensionalities, D, the accordance with 
HTSE results is better than that for the most critical for 
the SCGA S = 1/2 Ising model (n = D = 1). One 
can see from the Table I, that the SCGA yields for the 
0(n) models with n < 10 more accurate results than the 
1/n expansion performed-for the fee, bec and sc lattices 
by Okabe and MasutanCJ, who calculated numerically 
the analytical expressions (the double integcals-over the 
Brillouin zone) obtained by Abe and HikamieZrEj. 

The models with higher lattice dimensionalities, d, are 
very important for testing the present approximation, 
which allows for Gaussian, i.e., non-interacting, fluctu- 
ations in the system. The interaction between fluctua- 
tions dies-put in the spherical limit D — * oo, as well as 
for d ^ 4c3. Thus, one can expect that the SCGA yields 
rather accurate results for d ^ 4, whereas the deviations 
from the MFA described by the SCGA are still apprecia- 
ble. Indeed, from Fig. || one can see that for the 1(1/2) 
model on the d = 4 hypercubic lattice the distance be- 
tween upper and lower spinodal boundaries is very small. 
The latter is due to the fact t hat the singularity of the 
function P{X) at X — > 1 [see (4.2)] is only logarithmic. 
The discontinuity of magnetization at T c is for the 1(1/2) 
model still substantial, but decreases quickly with the in- 
crease of the number of spin components D. For d 5 
discontinuities of thermodynamic functions in the SCGA 
disappear. The value of 9 C for the 1(1/2) model in d — 4 
dimensions (see Table I) is over 2% iess then the 1/d ex- 
pansion result of Fisher and Gaunta, and this discrep- 
ancy diminishes smoothly with the increase of d. This 
can be seen from the comparison with the 9 C values of 
Ref. Hfor d — 5,6 and the recent high- accuracy results 
of Ref. jlfj for d = 6, 7. For d = 6, 7 the SCGA yields the 
6 C values 0.894 [0.90227 (Ref. ||), 0.90290 (Ref. ^] and 
0.913 [0.91922 (Ref. [h])], respectively. These values of 9 C 
are also very close to those for the spherical model (see 
Table I). With the increase of the spin dimensionality 
the accuracy of the SCGA also increases. In Table I the 
SCGA-results for 9 C for d — 4,5 are compared with those 
of the general- ro 1/d expansion by Gerber and Fisherc3 
terminated by the term <i~ 5 . Unfortunately, the termi- 
nated 1/d expansion becomes less accurate for larger val- 
ues of n and does not reproduce, unlike the SCGA, the 
exact results for the spherical model. It should be also 
noted that for the O(n) models with d ^ 4 the results for 
9 C approach with the increase of n those for the spherical 
model much faster than in three dimensions. This means 
that the coefficient in the 1/n expansion for 9 C in Refs. 
0, H, EH should be very small in high dimensions. 

The values of the energy of 3-dimcnsional spin systems 
on the upper spinodal boundary of the SCGA equations 
are also rather close to the series ones. The calculated 
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TABLE I. The values of reduced Curie temperatures 8 C = T c /T^ iFA calculated for different classical spin models on different 
lattices from the upper spinodal boundary of the SCGA equations. The results of other methods are placed below for comparison. 



Model \ Lattice fee 


bee 


sc 


Diamond fec-sc 


bec-sc 


hpc(4d) 


hpc(5d) 


(z = 12) 


(z=8) 


(z=6) 


(z=4) (* = 18) 


(z = U) 


(z = 8) 


(2 = 10) 



0(1) 

(Ising, S=l/2) 



0(2) 

(plane 

rotator) 

0(3) 

(Heisenberg 
S = oo) 

0(4) 



0(5) 
0(6) 

0(8) 

0(10) 



0(oo) 
(spherical) 

Ising, S — l 



Ising 

(n = l, D = 2) 

Ising, S^oo 
(n=l, D=3) 



x-y, S = oo 
(n = 2, L> = 3) 



0.808 
0.8161 





81620L2 
81628E 



0.797 n 
0.8033LD 



0.790 n 
0.794M-, 

0.7943N 



0.785 



0.781,- 

0.7981 



0.77 
0.7891 



90 



07780 



°- 771 n 

0.771H 
0.74368 



0.844p 
0.8510 n 
0.85246- 



0.845 



0.868, 
0.874) 
0.87682 
0.87698 



0.785 □ 
0.79385EJ 



0.773 r 
0.780212 
0.78019t 




0.761 r-| 

0.76306Eil 



0.757,-, 
0.771E3 

0.754-. 

0.762E3 n 

0.75295EI 



0.749L-. 

0.75lK 
0.74640C 

0.74' 
0.74 
0.7418 

0.71777 
0.826 



0.827 



0.853 



40 



0.744 
0.7517: 
0.75101 

0.729 p 
0.7343O-. 
0.7338SO 
0.7332H 

0.719r-| 

0.723E|Li 

0.7216E3-, 

0.72148Eil 

0.712 r-| 

0.7123flEil 
0.7123EJ 

0.707,-, 

0.728EI 

0.704^ 

0.717EL 

0.70009EI 

0.698,—. 

0.702EL 

0.6922lEl 

0.69- 
0.69> 
O.68681 

0.65946 



0.790 n 
0.79893B 



0.792 



0.822 
0.83195 



0.673 p 
0.6760E 



0.651 



0.637 



0.628 



0.854 E 
0.8609E 



0.847 



0.842 r 
0.8505t 



0.839 



0.824 n 
0.8307a 



0.813 



0.807 n 
0.8145B 



0.803 



0.812 n 
0.8340lH 



0.808 
0.829 



20 



0.864 n 
0.87694^ 



0.810 n 0.864 
0.8319EI 0.8762 



0.864 
0.874 



70 



0.808 n 0.864 
0.8275^ p 0.8737 
0.821oMy 



0.828 n 
0.8354LD 



0.808 p-, 0.768 p 
0.8156EI 0.7760EZI 



0.621 


0.837 


0.799 


0.807 
0.8262 


a 


0.864 n 
0.8730EI 


0.616 


0.835 


0.797 


0.807 
0.8253 





0.864 n 
0.8724EJ 


0.608 


0.832 


0.793 


0.807 
0.8240 





0.864 n 
0.8717EJ 


0.603 


0.830 


0.790 


0.806 
0.8240 





0.864 1-, 
0.8711a 


0.55776 


0.81397 


0.76656 


0.8066 
0.8150 


4sl 


0.86483-, 
0.8674EI 


0.727 


0.883 


0.857 


0.853 




0.896 


0.730 


0.883 


0.858 


0.854 




0.896 


0.767 


0.902 


0.879 


0.880 




0.916 


0.699 


0.871 


0.842 


0.843 




0.890 
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FIG. 4. Temperature dependences of magnetization in 
magnetic field m(H,T) for the classical Heisenberg model in 
the SCGA, compared with MC-simulations of Ref. 
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FIG. 5. Temperature dependences of magnetization and 
zero-field susceptibility of EuO in the SCGA, compared with 
the neutron scattering data of Ref. 



values of the normalized energy U(9 C ) = U(9 C )/U(0) of 
the 5 = 1/2 Ising model are 0.276 (0.25), 0.298 (0.27), 
0.365 (0.33) for the fee, bee and sc lattices, where the 
HTSE results are placed in brackets for comparison. For 
the classical Heisenberg model the normalized critical en- 
ergies are given by 0.237 (0.245), 0.261 (0.265), 0.315 
(0.325), respectively. For systems with higher lattice di- 
mensionalities, d, the energies are close to the ones for the 
spherical model, U(8 C ) = 1 — 9 C , especially for systems 
with many spin components. The normalized critical en- 
ergies in the SCGA of the 1(1/2), H(oo), and spherical 
models, respectively, are 0.225, 0.198, 0.1933 for d = 4 
and 0.143, 0.135, 0.1352 for d = 5. 

In a magnetic field, H ^ 0, the SCGA becomes more 
accurate, because the system is driven away from the crit- 
ical point (0, T c ), where the SCGA breaks down. The lat- 
ter leads to the disappearance of the fictitious first-order 
phase transition in the SCGA starting from the fields, 
which are much smaller than the exchange interaction 
(i.e., for h = H/Jo -C 1). For systems with a continuous 
spin symmetry (e.g., for the isotropic Heisenberg model) 
the magnetic field introduces a gap in the spin-wave spec- 
trum and suppresses the singular contribution to magne- 
tization oc (9 3 / 2 [see (4.7) and the following discussion], 
which improves the situation in the whole region below 
T c . A comparison of the SCGA results for the magnetiza- 
tion in magnetic field m(H, T) of the classical Heisenberg 
model on the sc lattice with tha^lC-simulation results 
of Binder and Muller-KrumbhaartiH is represented in Fig. 

I 

By application of the SCGA to experimentally inves- 
tigated magnetic systems one should resrict oneself to 
the ones with large spin values (5 3> 1) and to the 
temperature range T > T c /S, where the whole Bril- 
louin zone is populated by spin waves and the system 
behaves classically. An attempt to apply the SCGA to 
the 5=1/2 Heisenberg model using the Brillouin func- 
tion with 5 = 1/2 [i.e., the Langevin function ( |3.7[ with 
D = 1), which corresponds formally to the considera- 
tion of the model with n = 3 and D = 1, yields for the 



sc lattice, in addition to the wrong linear behavior of 
the magnetization at low temperatures, the phase tran- 
sition point 9 C = 0.592. being considerably higher than 
the HTSE value 0.56CO. On the other hand, for systems 
with 5^1 quantum effects in the range of elevated tem- 
peratures are determined by a small parameters 1 / (zS) 
and can be partially taken into account in the SCGA 
by using the Brillouin function Bg. In typical cases this 
introduces errors that are smaller than the intrinsic inac- 
curacy of the SCGA. The Heisenberg ferromagnets EuO 
(T c ~ 69 K) and EuS (T c ~ 16.6 K) having 5 = 7/2 are, 
perhaps, the most convenient materials for testing the. 
SCGA, and they were .extensively studied with NMRE3 
and neutron scatteringLJ methods. EuO and EuS form 
fee lattices, and the exchange interaction extends up to 
the next nearest sc neighbors. The contribution of dipole- 
dipole interaction (DDI) to T C MFA iM 1.7% for EuO and 
4.9% for EuS. With the use of HTSE it was shownEll that 
DDI suppreses to some extent the reduced transition tem- 
peratures 9 C = T C /T™ FA due to its competing nature. In 
the SCGA rigorous taking into account DDI requires al- 
lowing for correlations between different components of 
molecular field fluctuations, l a p with a ^ 13, which leads 
to the complication of the formalism and goes beyond 
the scope of this article. Instead of it, for a comparison 
with experimental data on EuO and EuS we take DDI 
into account in a simplified manner, only through the 
renormalization of T" c MFA mentioned above. The results 
of numerical calculations for EuO represented in Fig. ^ 
show the same accordance of the SCGA results with the 
neutron scattering data of Ref. [59] as its accordance with 
the HTSE and MC results demonstrated above. For EuS, 
due to the negative value of the n.n.n. exchange constant 
J2, and hence the reduction of the effective number of in- 
teracting neighbours, the level of fluctuations is greater 
than in EuO, the Watson integral W is close to that for 
the sc lattice, and the deviation of the results from the 
MFA , as well as the discrepancy between the SCGA and 
experimental results, is somewhat larger. 
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VI. DISCUSSION 

The self-consistent Gaussian approximation (SCGA) 
for classical spin systems described here is a unified the- 
ory applicable to a wide class of lattice models investi- 
gated currently by different groups with different meth- 
ods. The SCGA takes into account fluctuations of the 
molecular field in the simplest way and is sensitive to the 
lattice dimensionality and structure and to the form of 
spin interactions. The SCGA yields rather accurate val- 
ues of the field-dependent magnetization, m{H 1 T), and 
other thermodynamic functions in the whole plane (H, T) 
excluding the vicinity of the critical point (0, T c ). In par- 
ticular, the values of T c themselves can be determined 
in the SCGA with an accuracy better than 1%, which 
makes it already important for practical applications to 
new lattice and spin Hamiltonian types. 

Indeed, the SCGA is more flexible (although less ac- 
curate) than series expansions, and consideration of new 
substances reduces in the simplest case to som e modi- 
fications of the lattice integ ral P(X), Eq. (4.1), and of 
the Gaussian integrals ( 3.18 ). This can be exemplified by 
studying the crossover between fee, sc, and bec lattices 
varying the relative strength of the first and second near- 
est neighbor interactions Jij Ji, or the crossover between 
Ising, Heisenberg, and x-y models varying anisotropy 
const ants. Having solved the SCGA system of equations 
( |3.15 ), one obtains all quantities of interest as a result 
of a single calculation. More serious generalizations of 
the SCGA are required for consideration of systems with 
DDI or with a transverse field, where nondiagonal cor- 
relations of molecular field fluctuations should be taken 
into account. For systems with many interacting sub- 
lattices the number of variables in the SCGA equations 
increases quadratically with the number of sublattices, 
and calculations become cumbersome. 

Consideration of ferromagnets with a transverse field 
or antiferromagnets in field in the SCGA can be avoided, 
if one is interested only in zero-field susceptibilities. The 
zero-field ferro- and antiferromagnetic susceptibilities can 
be calculated through the correlation functions of the 
simplest ferromagnetic model with the longitudinal field. 
This requires, however, summation of some new diagram 
sequences and is the subject of a separate work. 

Possible improvements of the SCGA should include 
non-Ornstcin-Zcrnikc effects in spin correlation functions 
(CF's) and non-Gaussian fluctuations of the molecular 
field. The former seems to be more important, since using 
Ornstein-Zernike CF's leads, du e to singularities of the 
lattice integral P(X), Eq. ( |4.2| ), to the overestimation 
of fluctuational effects for 3-dimensional systems, which 
results in the breakdown of the SCGA in the critical re- 
gion. The diagram technique for classical spin systems 
used for the construction of the SCGA is undoubtedly 
the best instrument for its further development, because 
it allows summation of different more complicated dia- 
gram series than those considered here. All other pertur- 
bative schemes that do not take explicitly the advantage 
of classical properties of a system fail to reproduce the 



SCGA, although the physical picture of Gaussian fluctu- 
ations of the molecular field is quite transparent. 

One more possible application of the classical spin dia- 
gram technique is that to low-dimensional and finite-size 
systems, where the level of fluctuations is large and an 
improvement of the SCGA is necessary. The first step 
in this direction was the calculation of the energy and 
susceptibility of low-dimensipnal antiferromagnets in the 
whole temperature intervaEd and also for a nonzero mag- 
netic fieldE3 with the use of the 1/D expansion. By this 
calculation, the results of which are rather good even for 
D = 3, some diagram series going beyond the SCGA were 
summed up. This can, in principle, show how to improve 
the SCGA in a nonperturbative way with respect to D. 

The SCGA can be generalized also for inhomogeneous 
states of magnetics. It turns out, however, that interest- 
ing results can be o bta ined already in the limit D — > oo, 
where the model ([bl]) is analytically soluble but not 
equivalent to the standard spherical model of Berlin and 
Kadl3 in inhomogeneous situations, even in the is otro pic 



case. The anisotropic spherical model defined by (1.1) An. 
the limit D — > oo-was already applied to domain wallst3 
and to thin filmsE.1 

And, finally, of a principal importance would be to 
construct the dynamical part of the classical spin dia- 
gram technique and to try to generalize the SCGA for 
dynamics. 
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